1 Introduction

1.1 Research Questions

  1. Does the climate and vegetation satellite data make possible to get an accurate short-term prediction of agricultural drought prior the growing season onset over croplands areas in Chile?
  2. How variable is the short-term prediction accuracy at different time-lapse prior the growing season onset?
  3. How is the spatial variability of the variables importance for the model prediction across the croplands areas of Chile?

2 Methods

2.1 Unit under analysis

The census unit was selected mainly for the next reasons:

  1. The commune is the basic political-administrative unit for Chile and each commune is formed by several census unit.
  2. The census units are used for the agricultural census carry out each ten years, so there are several data about irrigation crops which allow characterizing each unit.
  3. Using census unit instead of commune allow to increased agricultural homogeneity and decrease the variability.
  4. This will help to the Ministry of Agriculture of Chile to choose a proper unit to be used for an emergency declaration by agricultural drought.

2.2 Summary

The aim of the proposal is to explore the short-term prediction of agricultural drought using remote sensing climate and vegetation data prior the growing season’s onset on croplands census units over Chile.

First, the delimitation of the study area using MODIS product is proposed. Landcover product MCD12Q1 allow us to create a cropland mask for Chile. The product MCD12Q2 of land cover dynamics will be used to calculate the onset of the greenup and dormancy for the growing season across the croplands areas.

For the prediction model, the output will be considered across the growing season and the input until six months before growing season.

The outcome proposed are:

  1. Mean standardized NDVI anomaly (\(\overline{zNDVI}_{gs}\))
  2. Minimum standardized NDVI anomaly (Maximum drought intensity, \(NDVI_{max_{gs}}\))
  3. Drought duration (\(DD_{gs}\))

The input for the prediction model will be considered from the onset of the growing seasons until six months before with a monthly periodicity. The variable proposed to be used in the model are:

  1. Standardized Precipitation Index (SPI) derived from CHIRPS 2.0 dataset and time-scales of 1, 3, 6, 12 months (5km).
  2. Monthly Evapotranspiration data (ET) from MOD16 product (1km).
  3. Soil moisture data from CCI products.
  4. Monthly Standardized NDVI anomaly (zNDVI).
Draft Sketch

Draft Sketch

Knowing that the last agricultural drought seasons over Chile were 2007-2008, 2008-2009 and 2014-2015; the following dataset are proposed to be used:

  1. 2007-2008 Training
  2. 2008-2009 Testing
  3. 2014-2015 Validation

2.3 Data

  1. Landcover product MCD12Q1 v5.1 IGBP scheme, 500m, yearly from 2001 to 2013.
  2. Land Cover Dynamic MCD12Q2 v5, 500m, yearly from 2001 to 2014.
  3. Vegetation Indices MOD13A1 v5 500m, 16 days, from 2000 to present.

3 Results

In this part the results obtained during the exploration of products MCD12Q1 and MCD12Q2 from sensor MODIS are presented.

3.1 Croplands Mask

3.1.1 Intensity of Croplands

To set the boundaries of the study area a cropland mask was created usign the layer which use the IGBP classification scheme from MCD12Q1 product. This products has yearly frequency and has data from 2001 to 2013. Considering the 13 years of data a intensity map was created which shows the 13 years percentege in which each pixel was classified as cropland (Fig. 1a). Map on Fig.1 compared information povided for the Agriculture Ministry of Chile about croplands (b) with that obtained from the IGBP scheme from MODIS.

Fig. 1: (a) Croplands intensity obtained from MCD12Q1 and (b) Croplands derived from information from the Chilean government.

Fig. 1: (a) Croplands intensity obtained from MCD12Q1 and (b) Croplands derived from information from the Chilean government.

3.1.2 Estimating a proper cropland mask

To estimate the intensity level to be used to derive a cropland mask, a validation was carried out to compare different cropland intensity levels against visual inspection of high-resolution imagery (google earth and ESRI imagery).

For the study area, 585 pixels samples (regular + random) were generated using as reference the pixels grid from the MCD12Q1 MODIS land cover product. Each sample has a dimension of 500m x 500m and was superposed over high-resolution images.

The IGBP classification scheme was used, this defines the cropland class as lands covered with temporary crops followed by harvest and a bare soil period, according to this definition the visual inspection was carried over each sample. For each sample, the proportion of cropland observed within its boundaries was estimated as a proportion between 0 and 100%.

In the cases where was not very clear the land cover type for the sample (mainly when could be confused with grassland) the inspection was made using historical google earth images, searching for bare soil periods as evidence of harvest.

The resulting validation data for the 585 samples is shown in the next dynamic map.

Then, were calculated ten cropland mask considering intensity from 10% to 100%, each 10%. Afterward, a confusion matrix was used for analyzing the intensity level and different statistics were derived . But, due that the samples data is imbalanced having 117 croplands and 468 no croplands, a Synthetic Minority Over-sampling Technique was used to transform the data to balanced data and then calculate the statistics, this was made over each intensity level (10%,20%, …,100%), and was repeated 10 times and averaged for each intensity level to have more stable results. This way it is possible to determine which intensity level has a better representation of the cropland area based on the visual inspection made over the 585 samples.

Next Figure show the statistics value variation (y-axis) from 10% to 100% cropland intensity (x-axis).

Fig 2: Variation of statistics values for different intensity levels between 10% and 100%

Fig 2: Variation of statistics values for different intensity levels between 10% and 100%

Sensitivity is a measure of the proportion of croplands samples correctly identified by the level of cropland intensity (mask). This value has its maximum at 10% and constantly decrease until 100%. Indicating, that increasing the cropland intensity, which decreases the spatial extension of the mask, causes that the correctly identified by the mask (true positive) decrease and increase the ‘no croplands’ samples that were identified as ‘no croplands’ by the mask. Similarly, the specificity which is a measure of the proportion of samples that were correctly identified as ‘no croplands’ by the mask, increase, due that decreasing the spatial extension of the mask, increase the number of samples that could be correctly identified as ‘no cropland’.

The balanced accuracy which is the mean of specificity and sensitivity present the higher value 0.78 at 30% cropland mask level and then constantly decrease until reach 100%. Similarly, global accuracy, kappa, and F1 shows the higher value at 30% with 0.78, 0.57 and 0.75 respectively and then decreases until reach 100%.

Finally, a cropland mask was derived using data provided for the Ministry of Agriculture of Chile (CONAF). This mask was analysed against the 585 samples. The statistics obtained were 0.8, 0.8, 0.61 and 0.77 for balanced accuracy, global accuracy, kappa and F1 respectively.

3.2 Phenological Analysis

3.2.1 Seasons

In this part, the product MCD12Q2 version 5 from MODIS was used for the phenology analysis. This product has 500m spatial resolution and yearly frequency from 2001 to 2014

Also, this product has eight layers from which in this case for the phenology analysis the following layers were used:

  1. Onset_Greenness_Increase: Greenup measure in days since January 1, 2000
  2. Onset_Greenness_Minimum: Dormancy measure in days since January 1, 2000

Each of the layers has two bands, the first band correspond with the first detected season and the second to the second season (if exists). Special care must to had on the using of this product, because is best aligned with northern hemisphere dynamics..

For seasonality analysis, the phenology layers were masked using the cropland mask (threshold = 30%) . Analyzing the layers for greenup and dormancy was possible to detect seasonality differences according to the greenup/dormancy onset timing. Fig. 3 present the seasonality map and three classes of seasons were identified: 1) uni-modality, 2) bi-modality, and 3) early uni-modality. The uni-modality season has a bigger extension representing a higher proportion of croplands over this part of Chile.

To explore the NDVI time-series through the two extensions 5 pixels were extracted on each season class from the MOD13A1 NDVI product. Each time-series was smoothed using ‘lowess’ and then the multi-annual average was calculated. The final time-series for seasonality are presented in Fig. 4. It is possible to observe that the ‘uni-modality’ seasonality has a greenup onset which could be identified around August and a dromancy onset close to April. In the case of the ‘early’ season is markedly different, the onset of greenup is on April/May and the onset dromancy at end of year.

To analyze the start of the season (SOS) and end of the season (EOS), only the pixels which were identified as ‘uni-modality’ and ‘bi-modality’ were used (not used early season). This was made with the goal of reducing the variability for the later analysis at census unit, also, these pixels may be representative of more homogeneous agriculture type. The months of the start of the growing season are mainly between July and August as showed in Fig. 3. In the other hand, the end of the season was identified from February to May.

Fig. 3: a) Seasonality classes according with greenup/dormancy onset timing, b) start of the season (SOS) and c) end of the season (EOS).

Fig. 3: a) Seasonality classes according with greenup/dormancy onset timing, b) start of the season (SOS) and c) end of the season (EOS).

Fig. 4: Extraction of NDVI time-series on 5 pixels for class.

Fig. 4: Extraction of NDVI time-series on 5 pixels for class.

Fig. 4: Extraction of NDVI time-series on 5 pixels for class.

Fig. 4: Extraction of NDVI time-series on 5 pixels for class.

3.2.2 Census unit level

The cropland mask was used to calculate the surface of cropland and number of pixels has each census unit. Next map present the census units over the study area.

Fig. 5: Map of spatial variation of (a) cropland proportion, and (b) number of pixels, in each censust unit.

Fig. 5: Map of spatial variation of (a) cropland proportion, and (b) number of pixels, in each censust unit. 

## Warning: Removed 115 rows containing missing values (geom_point).

Then the census units having a cropland surface starting on 30% were chosen. To have a measure of the seasonality variability on each unit the standard deviation was calculated for the SOS and EOS as presented on Fig. 6. The 75% of the units have less than 29 days standard deviatioon for the SOS and 27 days in the case of EOS.

Fig. 6: Standard Deaviation for SOS and EOS over the census units selected.

Fig. 6: Standard Deaviation for SOS and EOS over the census units selected.

Next map shows the Quantile 25% for SOS and 75% for EOS, as the possible threshold to be used as the onsets for the growing season.

Fig. 7: Quartile .25 for SOS and .75 for EOS

Fig. 7: Quartile .25 for SOS and .75 for EOS

3.3 Input Data

## Using coddis as id variables

## Using coddis as id variables

## Using coddis as id variables

3.4 Output Data

## Using coddis as id variables